

data = load('OscillatoryHyperShockData.mat');

expIndex = 14;
shockIndices = [0,0,0,0; 0,0,0,0; 0,0,0,0; 0,0,0,0; 0,0,0,0; 5069,5319,6022,6331; 3337,3522,4202,4381;
    3747,4024,4907,5137; 4026, 4331, 5490, 5688; 4830,5034,7204,7418; 351,640,865,1095;904,1196,1845,2195;
    812,1000,1711,1982; 1120,1369,2266,2567; 1431,1676,2904,3169; 5124,5474,6117,6541; 5721, 6168, 7216, 7691;
    3631,4015,4799,5085; 3689,3929,4997,5424;
    9286,9636,9700,10070; 1213,1583,2193,2647; 906,1238,1820,2167; 7213,7492,8407,8780; 7518,7853,9015,9356;
    7085,7473,8108,8578;4505,4812,5415,5739;4811,5193,6013,6338;5108,5361,6604,6882;3751,4130,7410,7673];

rates = cell2mat(data.GrowthRate(expIndex));
times = cell2mat(data.Time(expIndex));
magnitude = data.ShockPer(expIndex,1);
period = data.ShockPer(expIndex,2);

[row, col] = find(not(isnan(rates)));
times = times(col); rates = rates(col);


len = EffectiveLengthNoSTD(times, rates);


inds = shockIndices(expIndex, :);
ind1 = inds(1); ind2 = inds(2); ind3 = inds(3); ind4 = inds(4);

% strain1 = (l(ind1) - l(ind2))/l(ind1); strain2 = (l(ind4) - l(ind3))/l(ind3);
% ymodulus = strain1/strain2;
% interval = ( times(ind4) - times(ind1));
% 
% [interval, ymodulus]
% strain1

R0 =  0.00044;
unperL1 = 3.09678*exp((times - 59.5047)*R0);

%for the 200 shock with 150 seconds, index 24
% R0 =  0.00030;
% unperL2 = 3.633*exp((times - 381)*R0);


%for the 300 shock, 360 seconds
% R1 =  0.00013;
% unperL2 = 3.182*exp((times - 193.1)*R1);

%for the 300 shock, 120 seconds
% R1 =  -0.00001;
% unperL2 = 3.4*exp((times - 193.1)*R1);

%for the 300 shock, 120 seconds
R1 =  0.00002;
unperL2 = 3.35*exp((times - 193.1)*R1);

% R0 =  0.00034;%mean(rates(1:ind1))
% unperL = 3*exp((times)*R0);


strain1 = (unperL1(ind2) - len(ind2))/len(ind2); strain2 = (unperL1(ind3) - len(ind3))/len(ind3);
ymodulus = strain1/strain2;
interval = ( times(ind3) - times(ind2));

[R0*interval, ymodulus]
strain1


%%
figure('DefaultAxesFontSize',18)
% plot(times,l)

plot3(times,len,1:length(times), 'linewidth', 2); 

hold all

plot3(times,unperL1,1:length(times), 'linewidth', 2); 
plot3(times,unperL2,1:length(times), 'linewidth', 2);  
%plot3(times,rates,1:length(times)); 
view(0,90)

title(['magnitude: ' num2str(magnitude) ' , period: ' num2str(period)])

